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Abstract 

It is shown that the energy levels of the one-dimensional nonlinear Schrodinger, 
or Gross-Pitaevskii, equation with the homogeneous trap potential x 2p , p > 1, 
obey an approximate scaling law and as a consequence the energy increases ap- 
proximately linearly with the quantum number. Moreover, for a quadratic trap, 
p = 1, the rate of increase of energy with the quantum number is independent 
of the nonlinearity: this prediction is confirmed with numerical calculations. It 
is also shown that the energy levels computed using a variational approximation 
do not satisfy this scaling law. 

1 Introduction 

The Bose-Einstein condensate is described, approximately, by a mean-field approx- 
imation, see for example Friedrich (1998), that gives the Gross-Pitaevskii equation. 
In the one-dimensional problem considered here this equation takes the form, 

- r ^~ + \^ 2 x 2 y + A\y\ 2 y^Ey 1 (1) 

where x is the spatial coordinate, /x the atomic mass of the atoms comprising the 
condensate, u the classical frequency of a single atom in the trap potential. The 
nonlinear parameter A results from the use of a mean-field approximation to de- 
scribe the particle interactions and is defined in terms of fundamental constants, 
A = Airf^aaN/ \x where ao is the scattering length and N the effective density of 
atoms along the condensate axis. In most experimental circumstances the nonlinear 
constant A is large so perturbation methods are of little value. For the ground state, 
because the wave function varies relatively slowly and because the nonlinearity is 
large the Thomas-Fermi approximations, equation [n] below, provides a reasonable 
approximation to both the energy level and the wave function. For excited states no 
such simple approximation seems to be available. Yabulov et al (1997) have derived 



a re-normalised perturbation theory that gives approximate energy levels and wave 
functions, but we show in section || that this method seems to provide a poor estimate 
of the excited energy levels. 

In this paper we show that the energy levels satisfy a simple approximate scaling 
law and consequently that they are given approximately by the simple formula, 



E n (A) = 




The first term is just the Thomas-Fermi estimate of the ground state energy, obtained 
by neglecting the kinetic energy term. The second term is the dominant correction 
and is linear in n independent of A. We show also that the latter behaviour is a 
consequence of the particular form of the trap potential. 



2 Theory 

The eigenvalues of equation [j], E n (A), n = 0, 1, 2 • ■ •, are those values of E for which 
y(x) satisfy the boundary conditions \y\ — > as |x| — > oo and the normalisation 
condition 

p oo 

dx|y(ar)l 2 = 1. (3) 



For real eigenvalues we may assume y{x) to be real. 

Two of the four independent parameters in this equation may be removed by 
rescaling x and y and ensuring that the normalisation conditions is invariant, 

/ y' ^' * a i h 
i = ai, y = — = , u = — , A — aA , a = , 

rl y/Jl 

which replaces ji and fi by unity. In the following we drop all primes. 
By re- writing equation |l| in the form 



rftn, f)v 1 1 

dJ + W = ' v ^ = E ^y 2 Ay ' = 2 c 



and treating x as the 'time' we may interpret equation [j] as that of a classical particle 
of unit mass moving in a time-dependent potential, V(y, x). Conventional methods 
of classical dynamics provide a means of estimating the eigenvalues. 

The potential V(y, x) is stationary at y — and this is a minimum for times 
x < xq = y/2E/u> and for these times there are also maxima at 

y 2 = y m (xf = E(x)/A. 

For larger times, when E(x) < 0, there is only a maximum at y = 0. Hence quasi- 
periodic motion is possible for small times but for larger times almost all orbits 
diverge as |x| — > oo: for every E > 0, however, there are initial conditions for which 
y(x) — > as x — > oo. 

To be specific consider the even solution with initial conditions y(0) = a > and 
y'(0) = 0. For small a and large enough E this orbit will oscillate in the potential well 



2 



until the barrier at y — y m (x) is low enough for the orbit to either escape or to ride 
on the barrier top and eventually to zero: most orbits escape to infinity. Examples 
of these types of orbit are shown in the following figure. Here E = 15.0810, A = 100, 
uj = 1 and a — ai — 0.23975967 and a = a\ ± 0.0000001; the converged solution is 
not normalised. 




Figure 1 Some examples of even solutions of equation [lj with E = 15.0810 
and j/(0) = a, given in the text. 

This figure shows that the required solutions with y(x) — > as |ar| — * oo comprise 
a quasi-periodic part, for |x| < Xt where Xt is defined in equation ^| below, and 
a monotonically decreasing segment for \x\ > Xt- It also shows that the distance 
between nodes is almost constant: reasons for this are discussed later. 

Consider the oscillatory region. When E ^constant it follows from the definition 
of the Jacobi elliptic function that the odd and even solutions are, respectively 



y = asn(r,k), (y(0) = 0), y = am(K — T,k), (y(0) = a), 
where K = K(k) is the complete elliptic integral of the first kind and 



(5) 



r = xV2E- 
The period of these oscillations is 

T = 



j2 4 k 2 = ^= 



Aa 2 



2E - Aa 2 



\JlE - Aa 2 



K(k). 



(6) 



When E is constant the action of the above oscillatory solution may be written in 
the form 



I 



1 



2EF{k), F(k) 



3nk 2 Vl + k 2 



((l + k 2 )E(k)-(l-k 2 )K(k)), (7) 



where E(k) is the complete elliptic integral of the second kind, and is not to be 
confused with the energy. For each E there is bound motion if < Aa 2 < E and as 
k increases from zero to unity F(k) decreases from 1 to 4y2/(3n) ~ 0.6. The action 
is bounded by < I < I s , where I s is the action of the bound, non-periodic motion 
on the separatrix, where Aa 2 — E [k = 1), 



4 -3/^ 



3vrA 



(8) 
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Now consider the effect of E decreasing, but changing little during one period 
of the unperturbed motion. The principle of adiabatic invariance (Percival and 
Richards, 1982, chapter 9) shows that the action is almost invariant. The separatrix 
action, however, is not constant and decreases to zero at x = xq where u>xq — s/2E. 
All orbits cease to oscillate before this time and if the change in E is sufficiently slow 
this change occurs when the action equals the separatrix action. If xt is this time it 
is given by the solution of 



4 / 1 



3/2 



^[E--^l) =AI(E) (9) 

where the action is evaluated at E, the initial value of E. Adiabatic invariance 
shows that the solution oscillates with a local period, T, given by equation ^, which 
depends upon x. However, the period although singular at E(x) = Aa 2 , does not 
change significantly until E(x) is close to Aa 2 , so the nodes of the wave function are 
almost equally spaced. 

The quantum number, n, that labels the state is the number of zeros in the 
eigenfunction. The ground state, n = 0, has no zeros: the first excited state is odd 
and has one zero at the origin and the second excited state is even and has two zeros. 
Thus the oscillatory parts of the solution are represented by orbits that encircle the 
phase-space origin (n + l)/4 times before approaching the origin almost parallel to 
the y'-axis. There are n/4 oscillations in the interval < x < x t so we have the 
approximate relation xt = nT/4. For later use it is convenient to introduce the 
scaled variables 

tt E 2E 

N=-un, £ = -, and z = > 2, 

in terms of which k 2 — l/(z — 1) and the quantisation condition becomes 

wx t (£,z) = J^g{z), g(z) = 2 ^- r - (10) 

For large z, g(z) = 1 + ^ + 0{z- 2 ). 

Finally, we need an approximation to the motion for x > X\. The value of y(xt) 
must be close to the barrier height, y{xt) — y m {xt)'- if y{xt) <C y m (x t ) the orbit 
would complete another \ period and if y(x t ) > y m {xt) it would escape. But if 
y(xt) — y m (xt) the required subsequent orbit is approximated by expanding about 
the point in phase space that follows the potential maximum, by making the canonical 
transformation 

y = Q + ym{x), — = P + -r~ 
ax ax 

and expanding the equations of motion to second-order. Then if xq > is the time 
E(xq) = for xt < x < xq the equations of motion are 

d Q n dp ~t?i \^ d2 y™ d2 y™ Euj2 
— = p > — = E i x )Q - -nr " 



dx dx ' dx dx 2E{x)^AE{x) 

These equations may be solved numerically and it is seen that Q(x) remains small 
provided both |-P(:ct)| and \E{xt)Q{xt) — y'm(xt)\ are small or zero. As x — > xq the 
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solution diverges. However, over the interval of interest this expansion shows that an 
approximate solution is 



2/0) - Vm{x) = 



A ' 



x t < x < x a , E(xq) = 0. 



(11) 



This is, of course, the standard Thomas-Fermi approximation, obtained from equa- 
tion |l| by ignoring the kinetic energy term. 

Some idea of the accuracy of the approximations ^| and [ll] is given in the next 
figure comparing these with an exact solution. In this case E = 15, A — 100 which 
gives a = 0.23976 and x t = 2.7272. 




Figure 2 Graphs of the exact solution, the Thomas-Fermi solution [ll] for xt < x < xq and the 
adiabatic solution p| for < x < xt in the case E = 15 and A = 100, for which a = 0.23976 and 

x t = 2.7272. 

In the next section we use equations pL || and || to approximate the eigenvalues of 
equation [l] and to obtain an approximate scaling law. 



3 An approximate scaling law 

Here we show that the approximations described above may be used to derive an ap- 
proximate scaling law relating the energy, E, quantum number n and the nonlinearity 
parameter A by the single equation, 



£ = H (toAN- s/2 ) , E 



2E 



E 
N 



(12) 



for some function H . A consequence of this is that the energy levels behave like those 
of the linear oscillator in that the difference E n+ i(A) — E n (A) is almost indepenent 
of n and also of A. 

In order to derive this relation we first express z in terms of £ using the adiabatic 
and the quantisation conditions, equations ^| and [nj| respectively. These equations 
may be combined to give 



2%/2 

~3^r 



i - 



g(*) s 

4£ 2 



3/2 



F(k) 



i 



(13) 



which, in principle gives z{£). The behaviour of this function is shown in the next 
figure where 1/z is plotted as a function of £ = E/N . 



•5 



0.4 
0.35- 

0.3 
0.25- 

0.2 
0.15 



£ 



As 



0,g 



2 4 6 8 10 12 14 16 18 20 

Figure 3 Graph of l/z(£). 

1 and F — > 1, and so 2£ — > 1: in this limit, 
1 _ 2V2 (4£ 2 - I) 3 / 2 



3tt 



{2£f 



2£ ~ 1. 



As £ increases l/z(£) increases monotonically to 1/2. 

The normalisation condition, equation ||, can be written in the form 



rT/4 j-xo 

1 = 2na 2 / dxsn(r, k) 2 + 2 dx 
Jo Jx t 



E{x) 
A 



(14) 



The first of these integrals may be evaluated using relations given in Abramowitz and 
Stegun (1965, section 16.25), so we have 



2na 2 



\/2E - Aa 2 



K{k)-E{k) 
P 



ZujA 



(2E) 



3/2 _ x t 
3A 



(6E-u 



2„,2\ 



In terms of the scaled variables introduced in equation 10 this becomes 

3/4 K(k) - E{k) \ g{zf } 



3 At) 
2A 3 / 2 



= (2^) 3/2 1 + ^ ~ 



4£ \ir k 2 y/z{z - 1) 



16£ 3 



(15) 



(16) 



Since k 2 = l/(z — 1) and z is a function of £ through equation |l3], the right hand 
side of this equation depends only upon £ . Thus £ is a function only of the variable 
u;AN~ 3 / 2 , which is the scaling law |l2]. 

This analysis can be carried further with more approximations, but first we show 
the graph of the ratio 

R ( £ ) = j^fji ( 17 ) 

which is seen from equation [ll| and the fact that z — ► 2, tends to unity as £ — > 00. 
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Figure 4 Graphs of the ratio R(£), equation [l7| and the difference 100(i?(£ ) — Ri {£)). 

Expanding equation ^| in powers of 1/z gives 



27V3/2 



(2£) 



1 



16£ 3 



1+ 4z 



(18) 



An analysis of R(£) suggest that 1 — R(£ ) ~ for large £ , that in this range z ~ i 
and that z changes relatively slowly with £. Thus a simple approximation to this 
ratio is given by setting z equal to its asymptotic value, z = 2, to give 
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The graph of 100(i£(£) — R\(£)) is shown in figure [| and this demonstrates the 
accuracy of this simple approximation. 

On using R\ to approximate R(£) in equation |l6| and rearranging the equation 
we obtain 



E n (A) 



2/3 



7tt 



higher order terms. 



(19) 



1 /3Au)Y 

2 V^ - / ' 32" 

The first term in this equation is just the Thomas Fermi approximation, which follows 
from the normalisation condition, equation [l4|, by setting x t = 0. The second term 
increases linearly with n and, because the trap potential quadratic, is independent of 
A. Higher-order corrections come from the expansion about the asymptotic value of 
z and are complicated and not warranted because of other approximations made. 

The scaling law [l^ exists because the trap potential is homogeneous in x, so 
the adiabatic condition ^| may be expressed in terms of only two variables. For 
the quadratic potential these are £ = = 41 and z = and it is the form 
of these variables that gives the scaling law [l^ and ultimately the energy level [l9]. 
If the trap potential is (cux) 2p /2p, p > 1, the scaled energy may be taken to be 
£ = 2EN- 2 p/( 2 p+V and then the scaling law [l| becomes 

Alu 



E : 



2p 



and the energy levels become 
E n (A) 



1 



Alu 



2p 
2p + l 



2p+l 



1-KUJTL 



P 



Auj 



(20) 



When p = 1 this reduces to equation [ji 
upon the nonlinearity, A. 
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but when p ^ 1 the coefficient of n depends 
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4 Variational method 



Yukalov et al (1997) have used re- normalised perturbation theory to obtain analytic 
approximations to the energy levels of the 3<i nonlinear Schrodinger equation. Here 
we show that this method is equivalent to a Euler-Lagrange variational method and 
that the resulting energy levels of the excited states do not satisfy the scaling law 
described in equation Thus this method cannot be as accurate as implied by 
Yukalov et al (1997). 
With the Lagrangian 



1 fdy 



(21) 



and treating the energy as the Lagrange multiplier we see that the Euler-Lagrange 
equations with the functional and the constraint 



J[y] = / dx [L(y,y',x)-Ey 2 



dxy{x) 2 = 1, 



gives equation with /i = h = 1, and that the energy is then given by 



E = 



dx 



x 2 2 2 
— 10 X Z 

2 



Az 4 



(22) 



where z(x) is a solution of the Euler-Lagrange equation. For trial functions satisfying 
the normalisation condition we may use the simpler functional 



J[y] = / dxL(y,y',x). 



A natural trial function is 



z{x) = S j — H n (ax) exp 



1 2 2 

-—ax 
2 



h 2 n = 2 n n\^ 



where a is the variational parameter. Then the functional 23 becomes 



J(a) 



aA 
2^1 



I 711 -Ml 



dw iJ n ( , u;) 4 e 



(23) 



(24) 



(25) 



This is stationary so the appropriate value of a is given by the positive root of 
uj 2 AI n 



(2n+l)hl 







(a 2 ^ 


--) 




-I) 










a 2 1 



~)+Pn. (26) 



If A = these equations give the unperturbed energy levels and if A is small 
perturbation theory may be used to obtain the equivalent of Yukalov et al (1997), 
equation 44. For A > 1 and n — they give E = 0.677(wA) 2/3 which is 3.4% larger 
than the Thomas- Fermi energy, given by the first term in equation In this limit 
of large A perturbation theory may be used to give 



E n = -(2n+l)u 2 B 2 '* [l + JL + f- 



15 15 



B 



AI n 



1 



(2n+ \)oJ 2 hl ' 



(27) 
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It is also clear from equations that E/N depends only upon the variable z = 
AI n /((2n + l)/i^-y^)> which is different from the scaling law derived in the previous 
section. 

5 Numerical results 

In this section we compare the behaviour of the energy levels of equation [l], computed 
numerically, with the predictions of the above formula, equations [l^ and 

One method of numerically solving equation [j] is to perform a two-dimensional 
search in the (a, E) plane, where E is the energy and for even solutions y(0) — a > 
and for odd solutions y'(0) = a > 0. These solutions must a) satisfy the quantisation 
condition, b) tend to zero as x — ► oo and c) satisfy the normalisation condition. 
Since most solutions are unbounded this calculation is expedited by using a good 
first approximation, which is given by 



where y m (x) is the Thomas- Fermi solution defined in equation [ll] and f2 = 2ir/T 
where T is the period defined in equation ^. In practice the harmonic balance approx- 
imation il 2 = 2E—2a 2 A/2 was used for fi. The oscillatory part of this approximation 
has a slowly increasing amplitude in order that y{x) is continuous aX x — x t - 

This approximation has two free parameters, a and E, which were varied using the 
Marquardt algorithm to find values that simultaneously satisfied the normalisation 
condition || and the quantisation condition [l(]. For A = 200 this crude approximation 
gives a relative error of less than 1% for the ground state and 5% for the 16 th energy 
level. 

In the second stage of the calculation we use the energy E found above and vary 
a to find a value at which |y(xf)| < 5, for some small 6 and where X{ = 1.25a;o. This 
was achieved using a shooting algorithm that that varied a according to the value of 
y(xf). The solution obtained in this manner is not normalised, but we find that for 
small changes in E, J Q f dx y{x) 2 depends approximately linearly on E so it is possible 
to interpolate the energy to obtain values of (a,E) that give a correctly normalised 
solutions. 

In the following table are shown energy levels for A — 100 and 200. The exact 
numerical values are well approximated by the straight lines E n ~ 14.04 + 0.66n and 
E n ~ 22.40 + 0.74n, for A — 100 and 200 respectively, and the gradient of these 
lines is close to that predicted by equation |ll| The energy levels of the variational 
method do not behave in this manner, particularly for large A, and we conclude that 
the excited energy levels given by the re-normalised perturbation method used by 
Yukalov et al (1997) is not accurate for the one-dimensional nonlinear Schrodinger 
equation. 
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A = 100 



n 


2 4 6 


E n (numeric 
E n (equation 
E n (equation 


ill 

m 

2$ 


14.02 15.37 16.69 17.98 
14.12 15.47 16.84 18.20 
14.60 18.70 20.34 21.69 



A = 200 



n 


2 4 6 


E n (numeric 
E n (equation 
E n (equation 


m 

1) 


22.42 23.87 25.34 26.86 
22.41 23.77 25.13 26.49 
23.17 29.53 31.77 33.34 



6 Conclusions 

We have shown that the energy levels E n of the Gross-Pitaevskii equation [l] satisfy 
the approximate scaling law [1^, which relates the variables E, n, u), A in a single 
equation, which leads to the approximate energy levels |2| We have shown that other 
homogeneous trap potentials lead to similar scaling laws but only the energy levels 
of the quadratic trap have a coefficient of n that is independent of the nonlinear 
constant, see equation |2(]. It is also shown that the energy levels of the re-normalised 
perturbation method of Yukalov et al (1997) arc equivalent to a simple variational 
method and do not satisfy the scaling law derived here. 

The method used to derive these results involves interpreting the Gross-Pitaevskii 
equation as a mechanical system with a slowly varying potential, so that the idea 
of adiabatic invariance can be used. With this equivalence the spatial coordinate 
becomes the time, so the generalisation to the 2d- or 3d Gross-Pitaevskii equation is 
not apparent. For symmetric, many dimensional systems, however a similar approach 
may be possible though there are some problems with singularities at the origin that 
need to be resolved. 
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